#include "stdafx.h"
#include "ckRasterHistogram.h"
#include "ckRasterBand.h"
#include "ckRasterMask.h"
#include "math.h"

namespace ck
{
	namespace data
	{
		RasterHistogram::RasterHistogram()
		{
			m_nSize = 0;
			m_hist = NULL;
			m_histp = NULL;
			m_mean = 0.0;
			m_stdev = 0.0;
			m_stdev2 = 0.0;
		}

		RasterHistogram::RasterHistogram(z_uint32 nSize)
		{
			m_nSize = nSize;
			m_hist = new z_uint32[nSize];
			memset(m_hist, 0, nSize * sizeof(z_uint32));

			m_mean = 0.0;
			m_stdev = 0.0;
			m_stdev2 = 0.0;
		}

		RasterHistogram::~RasterHistogram()
		{
			if (m_hist)
			{
				delete[] m_hist;
				m_hist = NULL;
			}
			if (m_histp)
			{
				delete[] m_histp;
				m_histp = NULL;
			}
		}

		z_uint32 RasterHistogram::GetSize()
		{
			return m_nSize;
		}

		z_uint32 RasterHistogram::GetValue(z_uint32 i)
		{
			if (i >= m_nSize)
			{
				return 0;
			}

			return m_hist[i];
		}

		z_uint32* RasterHistogram::GetHist()
		{
			return m_hist;
		}

		double RasterHistogram::GetValuep(z_uint32 i)
		{
			if (i >= m_nSize)
			{
				return 0.0;
			}
			if (!m_histp)
			{
				return 0.0;
			}

			return m_histp[i];

		}

		double* RasterHistogram::GetHistp()
		{
			return m_histp;
		}

		RESULTCODE RasterHistogram::Normalize()
		{
			if (m_histp)
			{
				delete[] m_histp;
				m_histp = NULL;
			}

			m_histp = new double[m_nSize];
			memset(m_histp, 0, sizeof(double) * m_nSize);

			z_uint32 *ptr = m_hist;
			double *pfr = m_histp;
			double nCount = GetCount();
			for (z_uint32 i = 0; i < m_nSize; i++, ptr++, pfr++)
			{
				*pfr = *ptr / nCount;
			}

			return RC_SUCCESS;
		}

		RESULTCODE RasterHistogram::Create(RasterBand* pBand)
		{
			RESULTCODE rc = RC_SUCCESS;
			size_t nPixelSize = pBand->GetPixelSize();
			m_nGreys = pow(2, (nPixelSize << 3)) - 1;
			ckPixelType type = pBand->GetPixelType();
			switch (type)
			{
			case ckPixelTypeByte:
				rc = Create_Byte(pBand);
				break;
			case ckPixelTypeUInt16:
				rc = Create_UInt16(pBand);
				break;
			}

			return rc;
		}

		RESULTCODE RasterHistogram::Create_Byte(RasterBand* pBand)
		{
			RESULTCODE rc = RC_SUCCESS;

			RasterMask* pMask = pBand->GetMask();
			z_byte* data = (z_byte*)pBand->GetData();
			bool hasData = (data != NULL);

			if (!hasData)
			{
				pBand->ReadData();
			}

			ck_raster_stats stats;
			pBand->ComputeStatistics(&stats);
			m_nSize = (z_uint32)stats.maxv + 1;
			CreateHist();

			z_uint32 w = pBand->GetWidth();
			z_uint32 h = pBand->GetHeight();
			z_uint32 i, j;
			z_uint32 v = 0;

			z_byte* ptr = (z_byte*)pBand->GetData();
			for (i = 0; i < h; i++)
			{
				for (j = 0; j < w; j++, ptr++)
				{
					if (pMask)
					{
						if (pMask->IsValid(i, j))
						{
							v = *ptr;
							m_hist[v] += 1;
						}
					}
					else
					{
						v = *ptr;
						m_hist[v] += 1;
					}
				}
			}

			if (!hasData)
			{
				pBand->ReleaseData();
			}
			return rc;
		}

		RESULTCODE RasterHistogram::Create_UInt16(RasterBand* pBand)
		{
			RESULTCODE rc = RC_SUCCESS;

			RasterMask* pMask = pBand->GetMask();
			z_uint16* data = (z_uint16*)pBand->GetData();
			bool hasData = (data != NULL);

			if (!hasData)
			{
				pBand->ReadData();
			}

			ck_raster_stats stats;
			pBand->ComputeStatistics(&stats);
			m_nSize = (z_uint32)stats.maxv;

			z_uint32 w = pBand->GetWidth();
			z_uint32 h = pBand->GetHeight();
			z_uint32 i, j;
			z_uint32 v = 0;

			z_uint16* ptr = data;
			for (i = 0; i < h; i++)
			{
				for (j = 0; j < w; j++, ptr++)
				{
					if (pMask)
					{
						if (pMask->IsValid(i, j))
						{
							v = *ptr;
							m_hist[v] += 1;
						}
					}
					else
					{
						v = *ptr;
						m_hist[v] += 1;
					}
				}
			}

			if (!hasData)
			{
				pBand->ReleaseData();
			}
			return rc;
		}

		RESULTCODE RasterHistogram::CreateHist()
		{
			if (m_hist)
			{
				delete[] m_hist;
				m_hist = NULL;
			}
			m_hist = new z_uint32[m_nSize];
			memset(m_hist, 0, m_nSize * sizeof(z_uint32));
			return RC_SUCCESS;
		}

		z_uint32 RasterHistogram::GetCount()
		{
			if (!m_hist)
			{
				return 0;
			}

			z_uint32 nCount = 0;
			for (z_uint32 i = 0; i < m_nSize; i++)
			{
				nCount += m_hist[i];
			}
			return nCount;
		}

		RESULTCODE RasterHistogram::ComputeStatistics()
		{
			m_mean = ((double)GetCount()) / m_nGreys;

			double sum = 0.0;
			for (z_uint32 i = 0; i < m_nSize; i++)
			{
				sum += pow((m_hist[i] - m_mean), 2);
			}
			m_stdev2 = (sum / (m_nGreys - 1));
			m_stdev = sqrt(m_stdev2);

			return RC_SUCCESS;
		}

		double RasterHistogram::GetMean()
		{
			return m_mean;
		}

		double RasterHistogram::GetStdev()
		{
			return m_stdev;
		}

		double RasterHistogram::GetStdev2()
		{
			return m_stdev2;
		}

		double RasterHistogram::ComputeS()
		{
			double sum = 0.0;
			double d3 = pow(m_stdev, 3);

			double* pfr = m_histp;
			for (z_uint32 i = 0; i < m_nSize; i++, pfr++)
			{
				sum += pow((i - m_mean), 3)*(*pfr);
			}
			return sum / d3;
		}

		double RasterHistogram::ComputeK()
		{
			double sum = 0.0;
			double d4 = pow(m_stdev, 4);

			double* pfr = m_histp;
			for (z_uint32 i = 0; i < m_nSize; i++, pfr++)
			{
				sum += pow((i - m_mean), 4)*(*pfr);
			}
			return sum / d4;
		}

	}
}
